/*
 * To change this template, choose Tools | Templates
 * and open the template in the editor.
 */
package fvm;

import algebra.ElementVector;
import algebra.HMean;
import algebra.PMatrix_W0;
import algebra.Vector;
import java.awt.Dimension;
import java.awt.Point;
import java.util.ArrayList;
import kontisim.EulerTStep;
import kontisim.ODESystem;
import listener.ParentEventListener;

/**
 *
 * @author max
 */
public class FVModel implements ODESystem {

    FVModelData fvmd; //Modelldaten: FVElemente...
    ParentEventListener parenteventlistener; //Event Schnittstelle
    PMatrix_W0 mat; // Matrix
    ElementVector elemvec; // Vektor mit FVElementen
    Vector rvec; // Rechte Seite
    EulerTStep eulerstep;
    Point p;
    Dimension d;

    /*
     * Konstruktor
     */
    public FVModel() {


        this.p = new Point(0, 0);
        this.d = new Dimension();
    }

    //Funktioin zur Assemblierung der Matrix
    public void assemble() {

        //Eigener Thread damit GUI nicht lahmt
        Thread assemblethread = new Thread(new Runnable() {

            @Override
            public void run() {

                //Variablen initialisieren
                
                elemvec = new ElementVector(fvmd.getElements());
                mat = new PMatrix_W0(elemvec.size(), elemvec.size());
                rvec = new Vector(elemvec.size());

                //Statusbalken in GUI anzeigen
                
                parenteventlistener.Event_ProgBarRange(0, elemvec.size());
                parenteventlistener.Event_ProgBarValue(0);
                parenteventlistener.Event_ProgBarVisible(true);

                
                for (int i = 0; i < elemvec.size(); i++) {

                    // Zeile der Matrix bestimmen
                    
                    parenteventlistener.Event_ProgBarValue(i + 1);

                    FVElement middle = elemvec.get(i);

                    if (!middle.isboundary()) {

                        FVElement north = elemvec.get(i).getNorthNeigh();
                        FVElement south = elemvec.get(i).getSouthNeigh();
                        FVElement west = elemvec.get(i).getWestNeigh();
                        FVElement east = elemvec.get(i).getEastNeigh();

                        Double kxx_west = HMean.getHarmMean(west.getperm().getKxx(), middle.getperm().getKxx());
                        Double kxx_east = HMean.getHarmMean(east.getperm().getKxx(), middle.getperm().getKxx());

                        if (kxx_west.isNaN()) {
                            kxx_west = 0.;
                        }
                        if (kxx_east.isNaN()) {
                            kxx_east = 0.;
                        }

                        Double kxy_north = HMean.getHarmMean(north.getperm().getKxy(), middle.getperm().getKxy());
                        Double kxy_south = HMean.getHarmMean(south.getperm().getKxy(), middle.getperm().getKxy());

                        if (kxy_north.isNaN()) {
                            kxy_north = 0.;
                        }
                        if (kxy_south.isNaN()) {
                            kxy_south = 0.;
                        }

                        Double kyy_north = HMean.getHarmMean(north.getperm().getKyy(), middle.getperm().getKyy());
                        Double kyy_south = HMean.getHarmMean(south.getperm().getKyy(), middle.getperm().getKyy());

                        if (kyy_north.isNaN()) {
                            kyy_north = 0.;
                        }
                        if (kyy_south.isNaN()) {
                            kyy_south = 0.;
                        }

                        Double kyx_west = HMean.getHarmMean(west.getperm().getKxy(), middle.getperm().getKxy());
                        Double kyx_east = HMean.getHarmMean(east.getperm().getKxy(), middle.getperm().getKxy());

                        if (kyx_west.isNaN()) {
                            kyx_west = 0.;
                        }
                        if (kyx_east.isNaN()) {
                            kyx_east = 0.;
                        }

                        double dx_west = middle.getsupppoint().x - west.getsupppoint().x;
                        double dx_east = east.getsupppoint().x - middle.getsupppoint().x;

                        double dy_north = north.getsupppoint().y - middle.getsupppoint().y;
                        double dy_south = middle.getsupppoint().y - south.getsupppoint().y;

                        double diag = 0;

                        if (!west.isnoflow() && !west.isflow()) {
                            diag += kxx_west / dx_west;

                        }

                        if (!east.isnoflow() && !east.isflow()) {
                            diag += kxx_east / dx_east;

                        }

                        if (!north.isnoflow() && !north.isflow()) {
                            diag += kxy_north / dy_north;


                        }

                        if (!south.isnoflow() && !south.isflow()) {
                            diag += kxy_south / dy_south;

                        }
                        diag /= middle.getdx();

                        double diag2 = 0;

                        if (!west.isnoflow() && !west.isflow()) {
                            diag2 += kyx_west / dx_west;

                        }
                        if (!east.isnoflow() && !east.isflow()) {
                            diag2 += kyx_east / dx_east;
                        }
                        if (!north.isnoflow() && !north.isflow()) {
                            diag2 += kyy_north / dy_north;

                        }
                        if (!south.isnoflow() && !south.isflow()) {
                            diag2 += kyy_south / dy_south;
                        }
                        diag2 /= middle.getdy();

                        diag += diag2;

                        mat.set(i, i, diag);


                        if (north.isboundary()) {

                            if (north.isnoflow() || north.isflow()) {

                                if (north.isflow()) {

                                    rvec.set(elemvec.getElementNumber(middle), rvec.get(elemvec.getElementNumber(middle)) + north.gethk());

                                }

                            } else {
                                double tmp = north.gethk() * kyy_north / middle.getdy() / dy_north + north.gethk() * kxy_north / middle.getdx() / dy_north;

                                rvec.set(elemvec.getElementNumber(middle), rvec.get(elemvec.getElementNumber(middle)) + tmp);

                            }
                        } else {

                            double tmp = -kyy_north / middle.getdy() / dy_north - kxy_north / middle.getdx() / dy_north;
                            mat.set(i, elemvec.getElementNumber(north), tmp);


                        }

                        if (south.isboundary()) {

                            if (south.isnoflow() || south.isflow()) {

                                if (south.isflow()) {

                                    rvec.set(elemvec.getElementNumber(middle), rvec.get(elemvec.getElementNumber(middle)) + south.gethk());

                                }

                            } else {
                                double tmp = south.gethk() * kyy_north / middle.getdy() / dy_north + south.gethk() * kxy_north / middle.getdx() / dy_north;

                                rvec.set(elemvec.getElementNumber(middle), rvec.get(elemvec.getElementNumber(middle)) + tmp);
                            }
                        } else {

                            double tmp = -kyy_south / middle.getdy() / dy_south - kxy_south / middle.getdx() / dy_south;

                            mat.set(i, elemvec.getElementNumber(south), tmp);

                        }

                        if (east.isboundary()) {
                            if (east.isnoflow() || east.isflow()) {

                                if (east.isflow()) {

                                    rvec.set(elemvec.getElementNumber(middle), rvec.get(elemvec.getElementNumber(middle)) + east.gethk());



                                }


                            } else {
                                double tmp = east.gethk() * kxx_east / middle.getdx() / dx_east + east.gethk() * kyx_east / middle.getdy() / dx_east;

                                rvec.set(elemvec.getElementNumber(middle), rvec.get(elemvec.getElementNumber(middle)) + tmp);

                            }
                        } else {

                            double tmp = -kxx_east / middle.getdx() / dx_east - kyx_east / middle.getdy() / dx_east;

                            mat.set(i, elemvec.getElementNumber(east), tmp);
                        }

                        if (west.isboundary()) {

                            if (west.isnoflow() || west.isflow()) {


                                if (west.isflow()) {


                                    rvec.set(elemvec.getElementNumber(middle), rvec.get(elemvec.getElementNumber(middle)) + west.gethk());


                                }

                            } else {
                                double tmp = west.gethk() * kxx_west / middle.getdx() / dx_west + west.gethk() * kyx_west / middle.getdy() / dx_west;

                                rvec.set(elemvec.getElementNumber(middle), rvec.get(elemvec.getElementNumber(middle)) + tmp);

                            }
                        } else {

                            double tmp = -kxx_west / middle.getdx() / dx_west - kyx_west / middle.getdy() / dx_west;

                            mat.set(i, elemvec.getElementNumber(west), tmp);

                        }


                    }


                }

                //Fertig
                parenteventlistener.Event_Assembled();
                parenteventlistener.Event_ProgBarVisible(false);
            }
        });

      
        assemblethread.start();

    }

    /*
     * Simulation starten
     */
    public void StartSim(final double t, final double dt) {

        //Init
        
       

        for (int i = 0; i < elemvec.size(); i++) {

            //alle alten h's löschen
            elemvec.get(i).deltehs();
       if(!elemvec.get(i).isboundary())elemvec.get(i).seth0((fvmd.maxbc+fvmd.minbc)/2);
            

        }
        
        
        
                
        
        

        
        //Neuer Thread damit GUI nicht hängt
        Thread thread = new Thread(new Runnable() {

            public void run() {

                parenteventlistener.Event_getGUILocation();
                double lowerb = Double.POSITIVE_INFINITY;
                double upperb = Double.NEGATIVE_INFINITY;

                Vector vec = new Vector(elemvec.size());
                Vector vec_ne = new Vector(elemvec.size());
                Vector vec_S0 = new Vector(elemvec.size());


                for (int i = 0; i < vec.size(); i++) {

                    vec.set(i, elemvec.get(i).geth(0));
                    vec_ne.set(i, elemvec.get(i).getne());
                    vec_S0.set(i, elemvec.get(i).getS0());


                }

                eulerstep = new EulerTStep();


                int steps = (int) (t / dt);

                parenteventlistener.Event_ProgBarRange(0, steps);
                parenteventlistener.Event_ProgBarValue(0);
                parenteventlistener.Event_ProgBarVisible(true);

                for (int i = 0; i < steps; i++) {

                    parenteventlistener.Event_ProgBarValue(i + 1);

                    vec = eulerstep.TimeStep(getModel(), 0, dt, vec, vec_ne, vec_S0);

                    if (i == steps - 1) {
                        for (int j = 0; j < vec.size(); j++) {

                            elemvec.get(j).addh(vec.get(j));

                        }
                    }

                    //Fluss ausrechnen

                    if (i == steps - 1) {

                        for (int j = 0; j < elemvec.size(); j++) {

                            FVElement elem = elemvec.get(j);

                            FVElement FLOWX = null;
                            FVElement FLOWY = null;

                            double qy = 0, qx = 0;

                            double qy1 = 0, qy2 = 0, qx1 = 0, qx2 = 0;

                            if (!elem.isnoflow()) {


                                FVElement north = elem.getNorthNeigh(), south = elem.getSouthNeigh(),
                                        west = elem.getWestNeigh(), east = elem.getEastNeigh();


                                if (!north.isnoflow() && !south.isnoflow() && !west.isnoflow() && !east.isnoflow()) {

                                    qy1 = HMean.getHarmMean(elem.getperm().getKyy(), south.getperm().getKyy()) * elem.getdx() * (elem.gethk() - south.gethk()) / ((elem.getdy() / 2) + (south.getdy() / 2));
                                    qy2 = HMean.getHarmMean(elem.getperm().getKyy(), north.getperm().getKyy()) * elem.getdx() * (-elem.gethk() + north.gethk()) / ((elem.getdy() / 2) + (north.getdy() / 2));

                                    double k = HMean.getHarmMean(elem.getperm().getKxx(), west.getperm().getKxx());
                                    double A = elem.getdy() * 1;
                                    double h1 = elem.gethk();
                                    double h2 = west.gethk();
                                    double dx = Math.abs(elem.getsupppoint().x - west.getsupppoint().x);

                                    qx1 = k * (h1 - h2) * A / dx;

                                    k = HMean.getHarmMean(elem.getperm().getKxx(), east.getperm().getKxx());
                                    A = elem.getdy() * 1;
                                    h1 = elem.gethk();
                                    h2 = east.gethk();
                                    dx = Math.abs(elem.getsupppoint().x - east.getsupppoint().x);

                                    qx2 = k * (h2 - h1) * A / dx;

                                    qx = qx1 + qx2;
                                    qy = qy1 + qy2;

                                }

//                              
                                elem.setq(1, -qx/2, qy/2);
                            } else {
                                elem.setq(1, 0, 0);
                            }
                        }
                    }

                }

                for (int j = 0; j < elemvec.size(); j++) {



                    if (elemvec.get(j).gethk() < lowerb && !elemvec.get(j).isnoflow()) {
                        lowerb = vec.get(j);
                    }
                    if (elemvec.get(j).gethk() > upperb) {
                        upperb = vec.get(j);
                    }
                }


                //Fertig
                parenteventlistener.Event_ProgBarVisible(false);
                parenteventlistener.Event_Sim_finished(null, lowerb - 1e-8, upperb + 1e-8, steps);


            }
        });

        thread.start();



    }
//Nördliche Randbediungen setzen
    public void setnorthbound(ArrayList<Integer> typs, ArrayList<Double> vals) {

        if (fvmd.getElements().get(0).size() - 1 > typs.size()) {
            for (int i = 1; i < fvmd.getElements().get(0).size() - 1; i++) {
                int j = (int) (((double) (i - 1) / (double) (fvmd.getElements().get(0).size() - 2)) * typs.size());
                int type = typs.get(j);
                double val = vals.get(j);
                FVElement elem = fvmd.getElement(0, i);
                if (type == 1) {

                    elem.setnoflow(false);
                    elem.setboundary(true, val);
                    
                    if(val<fvmd.minbc)fvmd.minbc = val;
                    if(val>fvmd.maxbc)fvmd.maxbc = val;


                }
                if (type == 2) {

                    elem.setflow(true);
                    elem.seth0(val);


                }
                if (type == 3) {

                    elem.setnoflow(true);
                }

            }

        } else {

            for (int i = 1; i < fvmd.getElements().get(0).size() - 1; i++) {

                int type = typs.get(i - 1);
                double val = vals.get(i - 1);
                
                if(val<fvmd.minbc)fvmd.minbc = val;
                    if(val>fvmd.maxbc)fvmd.maxbc = val;
                
                FVElement elem = fvmd.getElement(0, i);

                if (type == 1) {

                    elem.setnoflow(false);
                    elem.setboundary(true, val);


                }
                if (type == 2) {

                    elem.setflow(true);
                    elem.seth0(val);


                }
                if (type == 3) {

                    elem.setnoflow(true);
                }

            }

        }
    }
//Südliche Randbediungen setzen
    public void setsouthbound(ArrayList<Integer> typs, ArrayList<Double> vals) {

        if (fvmd.getElements().get(0).size() - 1 > typs.size()) {
            for (int i = 1; i < fvmd.getElements().get(0).size() - 1; i++) {
                int j = (int) (((double) (i - 1) / (double) (fvmd.getElements().get(0).size() - 2)) * typs.size());
                int type = typs.get(j);
                double val = vals.get(j);
                
                if(val<fvmd.minbc)fvmd.minbc = val;
                    if(val>fvmd.maxbc)fvmd.maxbc = val;
                
                FVElement elem = fvmd.getElement(fvmd.getElements().size() - 1, i);
                if (type == 1) {

                    elem.setnoflow(false);
                    elem.setboundary(true, val);


                }
                if (type == 2) {

                    elem.setflow(true);
                    elem.seth0(val);


                }
                if (type == 3) {

                    elem.setnoflow(true);
                }

            }

        } else {

            for (int i = 1; i < fvmd.getElements().get(0).size() - 1; i++) {

                int type = typs.get(i - 1);
                double val = vals.get(i - 1);
                
                if(val<fvmd.minbc)fvmd.minbc = val;
                    if(val>fvmd.maxbc)fvmd.maxbc = val;
                
                FVElement elem = fvmd.getElement(fvmd.getElements().size() - 1, i);

                if (type == 1) {

                    elem.setnoflow(false);
                    elem.setboundary(true, val);


                }
                if (type == 2) {

                    elem.setflow(true);
                    elem.seth0(val);


                }
                if (type == 3) {

                    elem.setnoflow(true);
                }

            }

        }
    }
//Östliche Randbediungen setzen
    public void seteastbound(ArrayList<Integer> typs, ArrayList<Double> vals) {

        if (fvmd.getElements().size() - 1 > typs.size()) {
            for (int i = 1; i < fvmd.getElements().size() - 1; i++) {
                int j = (int) (((double) (i - 1) / (double) (fvmd.getElements().size() - 2)) * typs.size());
                int type = typs.get(j);
                double val = vals.get(j);
                
                if(val<fvmd.minbc)fvmd.minbc = val;
                    if(val>fvmd.maxbc)fvmd.maxbc = val;
                
                FVElement elem = fvmd.getElement(i, fvmd.getElements().get(0).size() - 1);
                if (type == 1) {

                    elem.setnoflow(false);
                    elem.setboundary(true, val);


                }
                if (type == 2) {

                    elem.setflow(true);
                    elem.seth0(val);


                }
                if (type == 3) {

                    elem.setnoflow(true);
                }

            }

        } else {

            for (int i = 1; i < fvmd.getElements().size() - 1; i++) {

                int type = typs.get(i - 1);
                double val = vals.get(i - 1);
                
                if(val<fvmd.minbc)fvmd.minbc = val;
                    if(val>fvmd.maxbc)fvmd.maxbc = val;
                
                FVElement elem = fvmd.getElement(i, fvmd.getElements().get(0).size() - 1);

                if (type == 1) {

                    elem.setnoflow(false);
                    elem.setboundary(true, val);


                }
                if (type == 2) {

                    elem.setflow(true);
                    elem.seth0(val);


                }
                if (type == 3) {

                    elem.setnoflow(true);
                }

            }

        }
    }
//Westliche Randbediungen setzen
    public void setwestbound(ArrayList<Integer> typs, ArrayList<Double> vals) {

        if (fvmd.getElements().size() - 1 > typs.size()) {
            for (int i = 1; i < fvmd.getElements().size() - 1; i++) {
                int j = (int) (((double) (i - 1) / (double) (fvmd.getElements().size() - 2)) * typs.size());
                int type = typs.get(j);
                double val = vals.get(j);
                
                if(val<fvmd.minbc)fvmd.minbc = val;
                    if(val>fvmd.maxbc)fvmd.maxbc = val;
                
                FVElement elem = fvmd.getElement(i, 0);
                if (type == 1) {

                    elem.setnoflow(false);
                    elem.setboundary(true, val);


                }
                if (type == 2) {

                    elem.setflow(true);
                    elem.seth0(val);


                }
                if (type == 3) {

                    elem.setnoflow(true);
                }

            }

        } else {

            for (int i = 1; i < fvmd.getElements().size() - 1; i++) {

                int type = typs.get(i - 1);
                double val = vals.get(i - 1);
                
                if(val<fvmd.minbc)fvmd.minbc = val;
                    if(val>fvmd.maxbc)fvmd.maxbc = val;
                
                FVElement elem = fvmd.getElement(i, 0);

                if (type == 1) {

                    elem.setnoflow(false);
                    elem.setboundary(true, val);


                }
                if (type == 2) {

                    elem.setflow(true);
                    elem.seth0(val);


                }
                if (type == 3) {

                    elem.setnoflow(true);
                }

            }

        }
    }

    public FVModel getModel() {
        return this;
    }

    public ElementVector getelemvec() {
        return elemvec;
    }

    public void setLocation(Point p, Dimension d) {
        this.p = p;
        this.d = d;
    }

    public void addModelData(FVModelData fvmd) {
        this.fvmd = fvmd;
    }

    public void addParentEventListener(ParentEventListener parent) {
        parenteventlistener = parent;
    }

    public Vector getValue(double t, Vector h) {





        //S0 fehlt noch
        //ne fehlt noch




        return rvec.scalar(h).minus(mat.vec_scal_mul(h));

    }

 
}
